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Abstract 

In this paper we introduce a binary search algorithm that efficiently finds initial 
maximum likelihood estimates for sequential experiments where a binary response 
is modeled by a continuous factor. The problem is motivated by switching mea- 
surements on superconducting Josephson junctions. In this quantum mechanical 
experiment, the current is the factor controlled by the experimenter and a binary 
response indicating the presence or the absence of a voltage response is measured. 
The prior knowledge on the model parameters is typically poor, which may cause 
the common approaches of initial estimation to fail. The binary search algorithm is 
designed to work reliably even when the prior information is very poor. The prop- 
erties of the algorithm are studied in simulations and an advantage over the initial 
estimation with equally spaced factor levels is demonstrated. We also study the 
cost-efficiency of the binary search algorithm and find the approximately optimal 
number of measurements per stage when there is a cost related to the number of 
stages in the experiment. 

Key words: optimal design, binary search, logistic regression, complementary 
log-log, quantum physics, switching measurement 



1 Introduction 



Consider an experiment where the probability of success (or failure) is a mono- 
tonic function of a factor controlled by the experimenter. The functional form 
of the response curve is known but the parameters defining its location and 
slope are unknown. Our goal is to estimate these unknown parameters in the 
situation where the prior knowledge on the parameters is very poor. Intu- 
itively it is clear that the factor levels should be chosen such a way that both 
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successes and failures are measured as responses - otherwise the experiment 
provides very little information on the parameters of the interest. However, if 
the location and slope of the response curve are unknown, the optimal factor 
levels are also unknown. 

The problem described above is encountered in quantum mechanical exper- 
iments called switching measurements [6,5,12]. In the experiment, a Joseph- 
son junction circuit is placed in a dilution refrigerator at temperature below 
4.2 kelvin. The height of a current pulse measured in nanoamperes is the fac- 
tor controlled by the experimenter and the absence or presence of a voltage 
pulse is measured as a response. For a certain narrow range of current pulses, 
the presence of the voltage response is a random variable depending on the 
height of the current pulse, which allows the Josephson junctions to be used 
as ultra-sensitive sensors. Only imprecise limits of the sensitive range can be 
given in advance: the current must be above zero and a rough upper limit is 
obtained by measuring the resistance of the circuit in room temperature. 

In the literature on the optimal design for binary response models, it is of- 
ten assumed that prior estimates are readily available or that they can be 
easily achieved from a pilot experiment. A natural but not always very effi- 
cient approach to initial estimation is to choose the initial factor levels that 
are equally spaced on some initial range. This approach has been employed 
in several theoretical and practical works, e.g. [13,14,3,2]. A more advanced 
approach based on the minimax principle was proposed and studied by Sitter 
[13] and King and Wong [9] 

In this paper we study the problem of initial design in detail and introduce 
a binary search algorithm that reliably finds initial maximum likelihood esti- 
mates (MLEs) even in the case where the prior information on the parameters 
is very poor. The binary search algorithm was already successfully applied 
to switching measurements [8]. We also study the cost-efficiency of the binary 
search algorithm and find the approximately optimal number of measurements 
per stage assuming that there is a cost related to the number of stages in the 
experiment. 

The rest of the paper is organized as follows. In Section 2 we review some 
commonly used binary response models and present the related theory of 
optimal design. The problem of the non-existence of MLEs with small and 
moderate sample sizes is studied. In Section 3, the binary search algorithm 
is introduced and its properties are studied. In Section 4, the binary search 
algorithm is compared to the initial estimation with equally spaced factor 
levels. In Section 5, we consider a cost model where there is a cost related to 
the number of measurements and a cost related to the number of stages and 
find the approximately optimal number of measurements in the binary search. 
In Section 6, we apply the results to switching measurements and discuss the 
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use of binary search in applications similar to the sport fishing example [13,9]. 
Section 7 concludes the paper. 



2 Binary response models 



There are three parametric models, the logit, the probit and the complemen- 
tary log-log (cloglog) model, that are frequently used to model the dependence 
between a binary response Y and a continuous factor x. All these models can 
be presented in the framework of generalized linear models [10] 

P(Y = 1) = E(Y) = F(ax + b), (1) 

where the response curve F is a cumulative distribution function (cdf) and the 
a and b are the parameters of the model to be estimated. The three response 
curves commonly used are: logistic distribution for the logit model 

, . exr>(ax + b) 

F(ax + b) = ^- J — , 2 

V ; 1 + exp(a:r + b) ' V ; 



normal distribution for the probit model 

F(ax + b) = <$>(ax + b) = ^ -L exp f -( ax + b ) 2 \ dx? (3) 

J-00 \/2ii \ 2 / 



and the Gompertz distribution for the cloglog model 

F(ax + &) = !- exp(- exp(aa; + b)). (4) 



The theory of optimal design provides results describing how the factor levels 
should be chosen in order to estimate parameters a and b optimally. Table 1 
adapted from [4] presents the (locally) D-optimal factor levels for the logit, the 
probit and the cloglog model. The apparent problem discussed in more detail 
by Minkin [11] and Sitter [13] is that the optimal factor levels are functions 
of the unknown parameters that we should estimate. It follows that we can 
design the optimal experiment only if we already know the parameters we 
want to estimate, which is an impractical requirement. The problem is often 
evaded assuming that we already have good initial estimates but less often it 
is discussed how the good initial estimates are obtained and what follows if 
initial information is very poor. 

Next we study the problem of the non-existence of MLEs. We consider maxi- 
mum likelihood estimation in small experiments where the sample size n < 200 
and assume that the true values of a and b are known so that the actual (lo- 
cally) D-optimal design can be used. MLEs cannot be estimated from the 
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Table 1 

D-optimal factor levels for the models (2), (3) and (4). The D-optimal factor levels 
x\ and x 2 are solved from equations ax\ + b = z\ and ax 2 + b = z 2 . Columns F{z\) 
and F(z 2 ) report the cdf levels related to the optimal factors, i.e., the probabilities 
of response 1. 

model D-optimal factors 

fi 4 F(4) F(4) 

logit -1.543 1.543 0.176 0.824 
probit -1.138 1.138 0.128 0.872 
cloglog -1.338 0.980 0.231 0.930 

measured data if it is possible to classify the responses to l's and 0's on the 
basis of the factor levels, e.g. response 1 is obtained if x > 100 and response 
is obtained if x < 100. More formally, the existence of MLEs for the model (1) 
requires that [1] 

max(xj|yj = 0) > mm(xi\yi = 1) and 

m&x(xi\yi = 1) > mm(xi\yi = 0), (5) 

where yi is the measured response for the factor level Xi. This condition is 
fulfilled, for instance, if there are two such factor levels that for each of them 
both 0's and l's are measured as responses. It follows that the probability 
for the non-existence of MLEs under the D-optimal design is equivalent to 
the probability that the only 0's or l's are measured as responses in either 
at z\ or at z 2 . According to the elementary laws of probability calculus, this 
probability can be expressed as 

P(No MLE) = q{ + p{ + q s 2 + p s 2 - q{q s 2 - q[p s 2 - p\q s 2 - p s lP s 2 , (6) 

where pi = F{z\) and^2 = F(z 2 ) are read from Table 1, q± = 1—pi, <?2 = 1 — P2, 
s = n/2 and n is the sample size, which is assumed to be an even number. 
The probabilities for the non-existence of MLEs in D-optimal designs as a 
function of sample size are presented in Table 2. It can be seen that surprisingly 
large sample sizes are needed in order to avoid the risk that MLEs cannot be 
calculated from the sample. The conclusion here is that even if we are using 
the D-optimal design, the non-existence of MLEs is a problem that cannot 
be ignored when the sample size is small. In practical situations where the 
D-optimal design is not known the problem may be even more serious. 
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Table 2 

Probabilities for the non-existence of MLEs when the D-optimal factor levels from 
Table 1 are used. 



sample size 


logit 


probit 


cloglog 


4 


0.91584142 


0.95047339 


0.95394709 


10 


0.61553275 


0.75551315 


0.77862340 


20 


0.26765833 


0.44578006 


0.52280966 


40 


0.04117204 


0.12633768 


0.23974289 


100 


0.00012482 


0.00217817 


0.02698092 


200 


0.00000001 


0.00000237 


0.00072786 



3 Initial estimation with binary search 

We start our introduction to the binary search algorithm considering types of 
the prior information. There is virtually always some prior information on the 
parameters available although this prior information may be very inaccurate 
and difficult to express in mathematical form. For instance, the following types 
of prior information may be encountered: 

• Point estimates of the parameters from the previous study/studies. 

• Range of the possible parameter values. 

• Prior distribution of the parameters. 

• Range of sensible factor levels. 

We concentrate to prior information given as a range of factor levels because 
it is the simplest type of prior information and the other types can be trans- 
formed to this type if needed. It is actually difficult to find a practical situation 
where no initial interval could be specified. Value zero is often a theoretical 
lower bound and some extreme high value can serve as the upper bound. In 
the other words, we can always specify the initial interval in such a way that it 
contains all sensible factor levels with high probability. If the prior information 
is very inaccurate, the initial interval can become very wide. 

The key idea of the binary search algorithm is to measure at such factor 
levels that the condition for the existence of MLEs (5) is fulfilled as quickly as 
possible. The motivation for this is that after the MLEs have been found, we 
can continue the experiment using the design that is optimal for the current 
estimates of the parameters. The proposed method is similar to the binary 
search algorithm that finds a root of a continuous function. Instead of a root 
we are looking for a factor level where both 0's and l's can be measured as 
responses. At every step, we measure the responses at the middle point of the 
current search interval: If only 0's are measured, the middle point is taken as 
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the new starting point of the interval. If only l's are measured, the middle 
point is taken as the new end point of the interval. The first part of the binary 
search ends when both O's and l's are measured at same point x. Another 
factor level needed to guarantee the existence of MLEs is then found in the 
neighborhood of x. The proposed algorithm can be presented as follows: 

1. Use the previous knowledge to construct such an interval [x m i n , a; max ] that 
we can be sure that -F(x max ; a, b) — F(x m i n ; a, b) is close to 1. Constructing 
this kind of interval is usually possible even if we have very little knowledge 
on the parameters. 

2. Use binary search to find such a point x from the interval [x m m, x max ] that 
both O's and l's are measured as responses. In binary search, we measure 
the responses at the middle point of the current interval. If only O's are 
measured, the middle point is taken as the new starting point of the interval. 
If only l's are measured, the middle point is taken as the new end point of 
the interval. Let x be found as the middle point of the interval [xj,x u ], i.e. 
x = (xi + x u )/2. 

3. Define |e| = (x u — Xi)/A. The sign of e is determined according to the 
measured response: sign(e) = sign(0.5 — y), where y is the average response 
for x. If 0.5 — y = the sign of e is chosen randomly. 

4. Measure first at point x + e. If both O's and l's are not measured as re- 
sponses, measure also at point x — e. 

5. If MLEs exist for the data measured so far, proceed to the maximum like- 
lihood estimation. Otherwise, divide e by two and return to Step 4. 

The progress of the algorithm is illustrated in Figure 1. Step 2 requires that 
we measure at least two times at each point. When the binary search in Step 
2 converges we have found one of the two points required for the maximum 
likelihood estimation. The other point that is needed is then found in the 
neighborhood of point x using again binary search. The procedure works re- 
liably and efficiently because at every step of binary search, the length of 
the current interval is divided by two. Consequently, the procedure converges 
exponentially fast regardless of the choice of the initial interval. 

After it has been found that condition (5) holds and the MLEs therefore exist, 
there are two alternative ways to use the data. We may calculate the MLEs 
from all data that has been cumulated in the binary search algorithm (MLE 
method I) or use only the endpoints of the initial interval and the two points 
found in Steps 2 and 4 of the algorithm (MLE method II). Neither of these 
approaches is fully satisfactory because the data originate from a procedure 
with data dependent stopping rules. However, the purpose of the initial design 
is only to produce some initial estimates for the primary experiment, which 
means that any sensible estimate based on real data is better than an arbitrary 
guess. 
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X 

Fig. 1. Binary search algorithm in action. The response curve of the cloglog model 
with a = 1 and b = is plotted. The initial interval is defined as [—5.7, 14.3]. The 
vertical lines present the measured factor levels and the Roman numerals indicate 
the order of the measurements. The numbers below the Roman numerals report 
the number of responses 1 per the total measurements; e.g. 3/5 means that three 
responses 1 were measured out of five measurements. The two points needed for 
the maximum likelihood estimation were found at measurements IV and VII. The 
obtained initial estimates of the parameters of the cloglog model and their standard 
errors are a = 1.83 (0.88) and b = 1.88 (1.29). 

The performance of the proposed binary search algorithm is studied in sim- 
ulations where the quality of the prior information varies from poor to very 
poor. We fix the parameters a = 1 and b = and specify the initial interval for 
factor levels as [h — d,h], where d > 2h is the length of the initial interval (a 
simulation parameter) and h is generated from Uniform(4, d — 4) distribution. 
The results can be generalized because arbitrary parameters a and b can be 
returned to the canonical case a = 1 and b = by a linear transformation. 
Note that the initial intervals are specified in such way that the middle point 
of the response curve always belongs to the interval but its relative location 
with respect to the endpoints of the interval is unknown and uniformly dis- 
tributed. The interval lengths d G {10,15,20,50,100,200,1000} are used in 
the simulations. The number of measurements per stage n& is also a simula- 
tion parameter and has 23 values {2, 3,4,..., 200, 500, 1000}. 5000 simulation 
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runs are performed for each combination of simulation parameters. For each 
simulation run we record 

• K, the number of stages needed for the convergence of the binary search, 

• 3k, the Fisher information matrix after the initial design, 

• a and b, the parameter estimates after the initial design. 

The validity of the initial interval is always first checked measuring the re- 
sponses at the ends of the initial interval. These measurements are counted as 
the two first stages and the actual binary search starts from the stage three. 

Tables 3, 4 and 5 report the simulation results for the cloglog model. It can 
be seen from Table 3 that more stages are needed when the initial interval 
widened but the number of the required additional stages is small. In fact, 
if the length of the initial interval is multiplied by two, one additional stage 
is required because at each stage of the binary search algorithm, the current 
interval is halved. It follows that if the length of the initial interval is multiplied 
by ten, log 2 (10) ~ 3.32 additional stages are required, which is in agreement 
with the differences between the observed K for d = 10, d = 100, and d = 1000 
in Table 3. 

The estimates using the MLE method I and the MLE method II are presented 
in Tables 4 and 5, respectively. The averages of estimates a and b indicate that 
a and b are systemically overestimated with the method I. The bias is large for 
small stage sizes (average a = 2.56 for n k = 3) and small for large stage sizes 
(average a = 1.04 for n k = 100) but it seems to be independent of the length 
of the initial interval d. For n k — 2, the average bias is very large indicating 
that some estimates were very large. The reason for the bias is that there is 
only limited number of outcomes when the initial estimation may stop. With 
method II, a is underestimated when n k = 2. On the basis of these results, it is 
clear that the stage size n k = 2 should not be used. With method II, parameter 
a is slightly overestimated for n k > 2 when d — 10 and underestimated when 
d = 100 and d = 1000. Parameter b is slightly underestimated. Similarly to 
method I, the bias decreases when n k increases. 

There are several ways to deal with the bias. The safest way is to choose large 
n k . If small n k is preferred, the results in Tables 4 and 5 could be employed 
in making a and b unbiased but this is not necessarily a good solution for 
method I because the distribution of a is strongly skewed, which can be seen 
comparing the mean and the 95th percentile. Simpler approach is to use the 
estimates without a bias correction and expect that the sequential design will 
quickly correct the bias. An overestimated a causes the D-optimal factor levels 
calculated according to Table 1 shift towards the center of the response curve. 
This is not particularly harmful because both 0's and l's can be measured as 
responses in the central area of the response curve. 
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The Fisher information after convergence of the binary search is almost in- 
dependent of the length of the initial interval. This can be explained by the 
fact that when the initial interval is wide, the measurement points at the first 
stages fall to the tails of the response curve and thus provide very little infor- 
mation in terms of the Fisher information. The number of stages that provide 
the information is therefore nearly constant. The same fact explains why the 
Fisher information is almost the same for method I and method II. 

Table 3 

Simulation results on the number of stages needed for the convergence using the 
binary search algorithm under the cloglog model. The numbers reported are the 
mean, the 5th percentile and the 95th percentile calculated from 5000 simulation 
runs for different lengths of the initial interval d and number of measurements per 
stage rifc. 





d 


= 10 




d 


= 100 




d 


= 1000 






mean 


5% 


95% 


mean 


5% 


95% 


mean 


5% 


95% 


2 


7.53 


5 


11 


10.93 


8 


14 


14.26 


12 


18 


3 


6.53 


4 


9 


9.87 


8 


13 


13.27 


11 


16 


4 


6.09 


4 


8 


9.52 


8 


12 


12.93 


11 


16 


5 


5.81 


4 


8 


9.25 


7 


12 


12.69 


11 


16 


10 


5.11 


4 


7 


8.72 


7 


11 


12.19 


10 


15 


100 


4.25 


4 


5 


7.83 


6 


10 


11.35 


10 


14 


1000 


4.09 


4 


5 


7.28 


6 


9 


10.91 


9 


14 



Full simulation results including results for the logit model and the probit 
model can be obtained from the author by a request. The results for the logit 
model and the probit model are otherwise similar to the results for the cloglog 
model but the estimates of b are unbiased because of the symmetry of the 
response curves. 



4 Comparison with equally spaced factors 

A commonly used practice for initial design is to choose factor levels that are 
equally spaced in the design space. Typically, the number of levels is selected 
on ad-hoc basis. More systematic approaches were proposed by Sitter [13] and 
King and Wong [9] who selected the initial design using the minimax prin- 
ciple. It was found by them that the number of the initial levels should be 
increased when the uncertainty about the model parameters increases. How- 
ever, when the initial range is wide and the number of initial levels is large, 
most of the levels will lie far from the D-optimal levels and thus provide very 
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Table 4 

Simulation results for the cloglog model using the binary search algorithm and 
MLE method I. The numbers reported are the mean, the 5th percentile and the 
95th percentile calculated from 5000 simulation runs for different lengths of the 
initial interval d and number of measurements per stage n^. 



a, estimated slope parameter (a = 1) 







d = 10 






d = 100 




d 


1 AAA 

= 1000 




nk 


mean 


5th 


95th 


mean 


5th 


95th 


mean 


5th 


95th 


2 


199.58 


0.50 


28.79 


221.84 


0.64 


27.20 


23.22 


0.58 


33.49 


3 


2.56 


0.60 


5.11 


21.03 


0.50 


7.17 


8.30 


0.64 


5.86 


4 


1.69 


0.57 


3.98 


1.64 


0.55 


3.64 


1.65 


0.58 


3.56 


5 


1.47 


0.63 


3.00 


1.46 


0.57 


2.85 


1.46 


0.57 


2.95 


10 


1.20 


0.67 


2.07 


1.19 


0.68 


2.03 


1.19 


0.61 


2.05 


100 


1.04 


0.85 


1.29 


1.02 


0.78 


1.27 


1.01 


0.79 


1.25 


1000 


1.00 


0.94 


1.07 


0.99 


0.83 


1.11 


1.00 


0.86 


1.15 


b, estimated location parameter 


(6 = 0) 
















d = 10 






d = 100 




d 


= 1000 




nk 


mean 


5th 


95th 


mean 


5th 


95th 


mean 


5th 


95th 


3 


0.67 


-1.05 


2.76 


40.77 


-1.05 


2.74 


14.72 


-1.01 


2.69 


4 


0.27 


-0.89 


1.78 


0.29 


-0.85 


1.83 


0.26 


-0.83 


1.76 


5 


0.20 


-0.76 


1.47 


0.20 


-0.75 


1.49 


0.19 


-0.76 


1.46 


10 


0.09 


-0.56 


0.92 


0.08 


-0.58 


0.86 


0.08 


-0.57 


0.86 


100 


0.02 


-0.21 


0.26 


0.01 


-0.26 


0.29 


-0.02 


-0.30 


0.25 


1000 


0.00 


-0.07 


0.07 


-0.03 


-0.26 


0.11 


0.01 


-0.15 


0.22 


D, the 


square 


root of the determinant of the Fisher information 










d= 10 






d= 100 




d 


= 1000 




n k 


mean 


5th 


95th 


mean 


5th 


95th 


mean 


5th 


95th 


2 


1.73 


0.08 


4.96 


1.66 


0.09 


3.75 


1.64 


0.07 


4.37 


3 


2.69 


0.51 


6.44 


2.66 


0.51 


6.79 


2.73 


0.55 


6.09 


4 


3.76 


1.08 


7.60 


3.61 


1.10 


7.60 


3.67 


0.94 


8.24 


5 


4.56 


1.30 


9.45 


4.46 


1.52 


9.13 


4.53 


1.51 


9.52 


10 


8.33 


3.25 


15.21 


8.35 


3.95 


14.76 


8.34 


3.80 


15.03 


100 


68.58 


30.98 


101.87 


55.56 


20.57 


89.39 


55.60 


24.47 


89.33 


1000 


679.14 


425.46 


952.25 


370.71 


103.53 


756.21 


323.80 


86.97 


658.73 



little information on the parameters of the interest. When the prior informa- 
tion is very poor, there is a substantial risk that equal spacing will just utilize 
resources without any guarantee of finding the MLEs. We study the perfor- 
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Table 5 

Simulation results for the cloglog model using the binary search algorithm and MLE 
method II. The numbers reported are the mean, the 5th percentile and the 95th 
percentile calculated from 5000 simulation runs for different lengths of the initial 
interval d and number of measurements per stage n^. 



a, estimated slope parameter (a = 1) 







d = 10 






d = 100 




d 


= 1000 




nk 


mean 


5th 


95th 


mean 


5th 


95th 


mean 


5th 


95th 


2 


0.86 


0.40 


1.70 


0.24 


0.09 


0.60 


0.06 


0.01 


0.18 


3 


1.21 


0.44 


3.19 


0.74 


0.10 


2.55 


0.78 


0.01 


2.04 


4 


1.04 


0.44 


2.52 


0.79 


0.11 


2.01 


0.73 


0.02 


1.80 


5 


1.02 


0.47 


2.26 


0.86 


0.11 


2.53 


0.80 


0.02 


2.02 


10 


1.03 


0.53 


1.87 


0.91 


0.23 


1.74 


0.90 


0.19 


1.73 


100 


1.02 


0.82 


1.28 


0.98 


0.71 


1.26 


0.96 


0.67 


1.24 


1000 


1.00 


0.94 


1.07 


0.98 


0.81 


1.11 


0.99 


0.84 


1.14 


b, estimated location parameter 


(6 = 0) 
















d = 10 






d = 100 




d 


= 1000 




nk 


mean 


5th 


95th 


mean 


5th 


95th 


mean 


5th 


95th 


2 


-0.02 


-0.96 


1.06 


-0.30 


-1.17 


0.36 


-0.39 


-1.23 


0.33 


3 


0.09 


-1.01 


1.20 


-0.13 


-1.47 


1.03 


-0.11 


-1.65 


1.06 


4 


0.01 


-0.92 


1.12 


-0.11 


-1.16 


1.02 


-0.15 


-1.23 


0.98 


5 


-0.00 


-0.91 


1.03 


-0.07 


-1.26 


1.01 


-0.12 


-1.45 


1.01 


10 


0.00 


-0.73 


0.80 


-0.08 


-1.00 


0.70 


-0.09 


-1.02 


0.71 


100 


-0.00 


-0.26 


0.26 


-0.02 


-0.37 


0.27 


-0.07 


-0.54 


0.24 


1000 


-0.00 


-0.07 


0.07 


-0.05 


-0.30 


0.10 


-0.00 


-0.22 


0.22 


D, the 


square 


root of the determinant of 


the Fisher information 










d = 10 






d = 100 




d 


= 1000 




nk 


mean 


5th 


95th 


mean 


5th 


95th 


mean 


5th 


95th 


2 


1.98 


0.22 


4.71 


3.28 


0.28 


9.77 


6.75 


0.34 


27.02 


3 


2.84 


0.42 


6.70 


3.87 


0.52 


12.32 


6.94 


0.65 


32.24 


4 


3.76 


0.91 


8.70 


4.20 


0.87 


12.92 


6.80 


0.78 


30.65 


5 


4.42 


1.08 


9.84 


4.42 


1.17 


13.16 


6.59 


1.14 


28.09 


10 


7.72 


2.96 


15.22 


6.81 


2.75 


13.10 


7.14 


2.76 


12.33 


100 


65.10 


27.30 


101.56 


50.68 


17.02 


83.85 


50.03 


17.46 


80.59 


1000 


670.40 


396.07 


952.25 


361.69 


87.16 


751.10 


314.30 


77.38 


652.07 



mance of initial estimation with equally spaced levels in the simulation setting 
presented in Section 3. The additional simulation parameters are the number 
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of levels M and the number of measurements per each level n m . M has values 
5, 10, 50, 100 and 1000 and n m has some suitably chosen values that lead to 
a reasonable sample size Mn m . For each combination of d, M and n m , 10000 
simulation runs are performed and the existence of MLEs is checked using the 
condition (5). 

Table 6 reports the empirical probabilities for the non-existence of MLEs under 
the cloglog model. The results show a dramatic change in the performance 
when the quality of the prior information decreases. For instance, if d = 10, 
the strategy with M — 10 and n m = 10 finds MLEs practically always but if 
d = 100 the same strategy fails almost without exceptions. If d = 1000, the 
existence of MLEs cannot be guaranteed even MM— 1000 and n m = 2. In 
general it seems that if the sample size Mn m is fixed, it is a better to have 
large M and small n m than vice versa. This is in agreement with the results 
in [13] and [9]. 

Initial estimation with equally spaced levels cannot compete with the binary 
search approach. Comparing Tables 3 and 6 we see that the required total 
sample sizes are significantly smaller in the binary search. By looking the row 
corresponding to n k = 5 in Table 3 we notice that using the binary search 
the sample size of 80 (= 5 • 16) was sufficient for finding the MLEs even when 
d = 1000. In Table 6, not even the sample size of 2000 (= 1000-2) was sufficient 
to eliminate the risk of MLE non-existence. 



5 Cost-efficient initial designs 

In this section we calculate cost-efficient stage sizes for the proposed binary 
search procedure. It is assumed that there is a cost related to the number 
of stages and a cost related to the number of observations. Without loss of 
generality, we fix the marginal cost of making one additional measurement to 
one and denote the cost of having one additional stage as C$- The same cost 
scheme was used for D-optimal sequential design in [7]. Our goal is to find the 
optimal number of measurements for initial estimation with the binary search 
algorithm as a function of the stage cost Cs and the length of the initial 
interval d. Optimality is defined here as the cost needed for finding the initial 
MLEs. We restrict only to initial designs were the number of measurements is 
the same for each stage of the binary search. If initial estimation with stage 
size n k takes K stages, the total cost becomes 

C(n k ) = KC s + Kn k . (7) 

When the stage cost Cs is fixed, we can find the stage size n k with the smallest 
total cost C{n k ) among the stage sizes used in the simulation in Section 3. 
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Table 6 

Empirical probabilities of MLE non-existence for the cloglog model in initial es- 
timation with equally spaced factor levels. The simulation parameters are d, the 
length of the initial interval, M, the number of factor levels, and n m , the number 
of measurements per each level. The results are means from 10000 simulation runs. 
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Using this approach the optimal stage size was calculated when the stage 
cost C s took 39 values from 0.001 to 1000. Based on these data, the following 
model was build for the the optimal rik 

log(rifc) = a + /31ogrf + 7logC 5 , (8) 

where the estimated parameters a, /3, and 7 for the logit, the probit and the 
cloglog model are presented in Table 7. In practical experiments, d is always 
unknown a priori and Cs may be known or unknown. We recommend that the 
choice of is based on the optimal rik calculated for several combinations of 
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d and Cs using equation (8). 
Table 7 

Summary of estimated models (8). 

logit R 2 = 0.9613 

estimate standard error 

Intercept a -0.964 0.069 

logd/0 -0.200 0.013 

logC S 7 0.839 0.012 



probit R 2 = 0.9765 

estimate standard error 
Intercept a -1.27 0.058 
logd/0 -0.223 0.011 

logC S 7 0.892 0.010 

cloglog R 2 = 0.9753 

estimate standard error 
Intercept a -1.063 0.057 
\ogd (3 -0.214 0.010 

logC 5 7 0.850 0.010 



6 Applications 



We reconsider the measurement data from [8] in order to illustrate the applica- 
tion of the binary search algorithm. In the experiment, a sample consisting of 
aluminium-aluminium oxide-aluminium Josephson junction circuit in a dilu- 
tion refrigerator at 20 millikelvin temperature was connected to computer con- 
trolled measurement electronics in order to apply the current pulses and record 
the resulting voltage pulses. The resistance of the sample at room temperature 
suggested that a pulse of 300 nA always causes a switching (response 1), which 
gave the upper limit for the initial estimation. The lower limit for the initial 
estimation, 200 nA was roughly estimated from the dimensions of the Joseph- 
son junction by an experienced physicist. Alternatively, nA may be used as 
the lower limit for the initial estimation. The initial estimation used the stage 
size 50. A major difference compared to the algorithm presented in Section 3 
was that the starting value of e was set as one. In this particular situation, 
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this was a good choice because the MLEs were found after measuring only at 
four points. Taking the final estimates a = 0.240 and b = —60.628 calculated 
from 117288 data points as the true parameters, we study the experiment in 
the light of the simulation results of Sections 3 and 5. The standardized length 
of the initial interval is d = 0.240 • (300 - 200) = 24, or d = 0.240 • 300 = 72 
for the wider interval. The standardized stage cost Cs was estimated from 
the data and the value Cs = 228 A was obtained [7]. The model (8) gives the 
optimal n k = exp(-1.063 - 0.214 log(24) + 0.85 log(228.4)) « 18 for the initial 
interval [200,300], and n k ss 14 for the initial interval [0,300]. According to 
Table 3 we expect that the number of stages needed is between 4 and 11. 

In addition to switching measurement, the binary search algorithm may be 
useful also in some other problems. Sitter [13] and King and Wong [9] con- 
sider an application where the problem was to evaluate the economic value 
of the sport fishing in British Columbia tidal waters. Fishermen were asked 
the question "If the cost of your fishing trip had been x dollars higher today, 
would you still have gone fishing?" Originally, x varied from $1 to $50 but 
in the middle of the survey it was found that a number of fishermen were 
willing to pay even more than $50 extra and the range of x was widened to 
cover values from $1 to $100. Obviously, the original survey was not designed 
optimally. As an improvement, Sitter [13] proposed a robust seven point de- 
sign with equal weights for this problem and King and Wong [9] improved 
Sitter's results by proposing a robust nine point design with unequal weights. 
This and similar surveys could benefit from sequential designs with efficient 
initial design. By using a cellular phone the interviewer could have a direct 
contact to the computer server collecting the data, which would allow the use 
of the binary search algorithm described in this paper. In this scenario, the 
interviewer sends the response of a fisherman to the server that after short 
calculations returns the amount of dollars for the next question. This would 
minimize the loss of efficiency due to badly chosen design space. 



7 Conclusion 

When the problem of optimal design in binary response experiments is consid- 
ered in statistical literature it is often assumed that the initial estimates are 
available. In this paper we propose a binary search algorithm that quickly and 
reliably returns initial estimates. After finding the initial MLEs the experiment 
may continue with likelihood based sequential estimation. 

Results provided in the paper consider the logit, the probit and the cloglog 
model but the proposed approach is not restricted to any particular paramet- 
ric model. On the other hand, the proposed approach relies on the idea of 
sequential design. If there is a considerable delay between choosing the factor 
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level and measuring the response, sequential designs may be impractical. This 
applies to many biomedical dose-response trials. In contrast, in many physi- 
cal experiments the response is measured almost instantly but the total time 
available for the experiment may be restricted. Efficient sequential designs 
may be very useful in these problems. 

The presented simulations compared the binary search algorithm to the initial 
estimation with equally spaced factor levels. The key difference between the 
approaches is that with equally spaced levels, a good guess for the initial 
interval must be available whereas with the binary search algorithm, any initial 
interval that contains the middle point of the response curve is sufficient. If 
the prior information is very poor, the proposed binary search approach has 
a huge advantage over the initial estimation with equally spaced levels. The 
importance of reliable initial estimation is emphasized by the fact that if the 
prior information is poor it is often difficult to deduce how poor it exactly is. 
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